function [P_apes, f] = apes(x, fs, M, NFFT)
% APES线谱估计函数
% 输入参数：
%   x     : 输入信号（一维向量）
%   fs    : 采样频率 (Hz)
%   M     : 滤波器长度（建议为信号长度的1/3到1/2）
%   NFFT  : 频率点数（默认1024）
% 输出参数：
%   P_apes: APES幅度谱估计值
%   f     : 对应的频率向量 (Hz)

if nargin < 4
    NFFT = 1024; % 默认频率点数
end

N = length(x);

% 1. 构造数据矩阵 H (Hankel矩阵)
L = N - M + 1; % 数据矩阵的列数
H = zeros(M, L);
for i = 1:L
    H(:, i) = x(i:i+M-1);
end

% 2. 计算数据协方差矩阵 R
R = (H * H') / L; % 协方差矩阵

% 添加正则化项确保矩阵可逆（可选）
epsilon = 1e-6 * trace(R)/M;
R = R + epsilon * eye(M);

% 3. 计算频率网格
f = linspace(-fs/2, fs/2-1, NFFT);
w = 2*pi*f/fs; % 归一化角频率

% 4. 计算APES谱
P_apes = zeros(NFFT, 1);
for k = 1:NFFT
    % 构造导向矢量
    a = exp(1j * w(k) * (0:M-1)');
    % 计算最优滤波器系数
    h = (R \ a) / (a' * (R \ a)); % h = R^{-1}a / (a^H R^{-1} a) -> w_{apes}
    % 计算APES幅度
    P_apes(k) = abs(h' * H * ones(L, 1)/L);
end

% 转换为dB尺度（可选）
% P_apes = 20*log10(abs(P_apes));
end